Introduction
Nous souhaitons réalisé l’étude d’une série temporelle et faire des prévisions sur celle-ci.
Cette série temporelle est le trafic mensuel d’une Compagnie aérienne de janvier 2011 à août 2019.
Nos prévisions portent sur les 8 mois de l’année 2019
Représentation graphique de la série.
Import des données
Import de la base, on sélectionne la colonne des valeurs de trafic.
library(readr)
data <- read_delim("Trafic-voyageurs.csv",
delim = ";", locale = locale(encoding = "ISO-8859-1"))data_value <- data[,2]
summary(data)## dates trafic
## Length:104 Min. :220876
## Class :character 1st Qu.:297154
## Mode :character Median :355178
## Mean :354651
## 3rd Qu.:407331
## Max. :505190Création de notre série temporelle
Nous représentons nos données sous forme de série temporelle.
Une série temporelle est un ensemble de métrique mesurée sur des intervalles de temps réguliers.
Création de la série chronologique avec la librairie TSstudio :
- start : année de début,
- frequency : La fréquence est mensuelle donc on sélectionne une fréquence de 12/an
library(TSstudio)
data_ts <- ts(data_value, start = 2011, frequency = 12)
plot_1_TimeSeries(data_ts)##
## Attachement du package : 'plotly'## L'objet suivant est masqué depuis 'package:ggplot2':
##
## last_plot## L'objet suivant est masqué depuis 'package:stats':
##
## filter## L'objet suivant est masqué depuis 'package:graphics':
##
## layoutSéparation jeux de données
data_ts_train <-
window(data_ts, start = c(2011, 1), end = c(2018, 12))
data_ts_test <- window(data_ts, start = c(2019, 1), end = c(2019, 8))
names(data)[1] <- "ds"
names(data)[2] <- "y"
data_train <- data[1:96, ]
data_test <- data[97:104, ]
plot(data_ts, xlim = c(2011, 2020))
lines(data_ts_test, col = 3)
legend(
"topleft",
lty = 1,
col = c(1, 3),
legend = c("Série chronologique Train", "Série chronologique Test")
)
On observe une forte tendance, et un paterne semble se répété, s’agit-il d’une saisonnalité ?
Analyse de notre série temporelle
Pour commencer nous allons nous intéresser à savoir si notre série est stationnaire ou non.
La stationnarité d’une série signifie que le processus qui génère la série ne change pas dans le temps. Cela ne veut pas dire que la série ne change pas dans le temps, mais que la façon dont elle change, n’est pas modifié dans le temps.
Testons si la série est stationnaire :
Par la suite la library tsutils permet de vérifier automatiquement s’il y a une tendance et une saisonnalité grâce au test effectué par la fonction seasplot().
library(tseries)
adf.test(data_ts) #p-value <0.5 => on ne rejete pas H0 => non stationnaire## Warning in adf.test(data_ts): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: data_ts
## Dickey-Fuller = -5.4857, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationarykpss.test(data_ts)## Warning in kpss.test(data_ts): p-value smaller than printed p-value##
## KPSS Test for Level Stationarity
##
## data: data_ts
## KPSS Level = 2.0803, Truncation lag parameter = 4, p-value = 0.01library(tsutils)## Registered S3 method overwritten by 'greybox':
## method from
## print.pcor lavaseasplot(data_ts) #TRUE | TRUE
## Results of statistical testing
## Evidence of trend: TRUE (pval: 0)
## Evidence of seasonality: TRUE (pval: 0)Donc notre série est bien non-stationnaire.
Décomposition de la série chronologique :
Généralement les séries chronologiques sont divisées en plusieurs composantes qui peuvent être additionnées ou multipliées entre elles.
Additive : Ts = Tendance + Saisonnalité+ Erreur
Multiplicative : Ts = Tendance * Saisonnalité * Erreur
Donc chaque point de notre série temporelle peut être exprimer par :
Tendance \(Tt\) : qui correspond à l’indicateur de si nos données changent sur le long terme ou non, elle permet de savoir si la série se déplace de manière fluide ou progressivement.
Saisonnalité \(St\) : décrit un pattern qui se reproduit sur un intervalle de temps régulier (ici une année), c’est à dire que chaque année, le même pattern est présent.
Erreur \(ϵt\) : ou résidus correspond aux “restes” quand on a supprimé la tendance et la saisonnalité. Ces valeurs peuvent s’expliquer par des évènements aléatoire qui affecte la valeur de la série (pandémie, choc boursiers, inflation …)
Pour la suite, nous supposons que la décomposition est additive.
decomposed_data <- decompose(data_ts_train, type="additive")
plot(decomposed_data$trend)
plot(decomposed_data$seasonal)
plot(decomposed_data$random)
boxplot(data_ts ~ cycle(data_ts))
Le choix d’une décomposition additive semble visuellement correct.
Analyse des résidus
L’ACF est le test de l’auto-corrélation par le décalage d’elle même, c’est un test important car il montre si les états précédents (observations décalées) de la série chronologique ont une influence sur l’état actuel.
Si l’auto-corrélation croise la ligne bleue, cela implique qu’un décalage est significativement corrélé avec la série actuelle.
On peut également regarder la distribution des résidus qui ici suit une courbe gaussienne centrée en 0 avec quelques outliers.
checkresiduals(remainder(decomposed_data))## Warning in modeldf.default(object): Could not find appropriate degrees of
## freedom for this model.
Représentation de la saisonnalité
En superposant les différentes années on peut observer une saisonnalité qui se répète chaque année.
Ce phénomène est encore plus remarquable quand on supprime la tendance.
On peut noter plusieurs “pic haut” :
Décembre / Janvier
Mars
Juin
Septembre (le plus haut)
Et plusieurs “pics faible” :
Avril
Juillet - Août (la période estivale)
Novembre
ggseasonplot(data_ts)
#on supprime la tendance
data_ts_without_trend = diff(data_ts)
SeasonPlot <- ggseasonplot(data_ts_without_trend) +
labs(
title = "Trafic sans la tendance",
subtitle = "Visualisation de la saisonnalité",
x = "Mois",
y = "Nombre de Voyageurs"
) +
geom_line(size = 1.1, alpha = 0.65) +
theme_fivethirtyeight() +
theme(axis.title = element_text()) +
scale_color_brewer(palette = "Paired") +
theme(axis.title = element_text(), text = element_text(family = "Rubik"))
SeasonPlot
Modèles Espace-état
Pour commencer nos prédictions des 8 premiers mois de l’année 2019, nous allons essayer quelques modèles basiques :
- meanf : Average Method : prend la valeur moyenne de toute les observations pour toutes les prédictions,
- naive : Naive Method : prend la dernière observation pour toutes les prédictions,
- drift : Drift Method : prend la première et la dernière observations et trace une lignes entre les deux, on utilise la courbe pour les prédictions,
- snaive : Seasonal Naive Forecast : Prend la dernière valeur de la saison précédente comme prédiction (ex : sept 2018 = sep 2019 + erreur)
Pour évaluer notre modèle on regarde :
MAE : Mean Absolute Error
RMSE : Root Mean Squarred Error
MASE : Mean Absolute Scaled Error
MAPE : Mean Absolute Percentage Error
Au formules associés :
res = pred - val
MAE = sum(abs(res))/length(val)
RSS = sum(res^2)
MSE = RSS/length(val)
RMSE = sqrt(MSE)
La plus populaire est la MAPE
MAPE(y_pred, y_true)
Concrètement si on a une MAPE de valeur 6% cela signifie la différence moyenne entre la prédiction et la valeur réelle est de 6%.
library(forecast)
mean <- meanf(data_ts_train, h=8)
naivem <- naive(data_ts_train, h=8)
driftm <- rwf(data_ts_train, h=8, drif=T)
snaivem <- snaive(data_ts_train, h=8)## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique## Warning in box(...): "plot.conf" n'est pas un paramètre graphique
## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique## Warning in box(...): "plot.conf" n'est pas un paramètre graphique
## Warning in plot.window(xlim, ylim, log, ...): "plot.conf" n'est pas un paramètre
## graphique## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "plot.conf" n'est
## pas un paramètre graphique## Warning in axis(1, ...): "plot.conf" n'est pas un paramètre graphique## Warning in axis(2, ...): "plot.conf" n'est pas un paramètre graphique## Warning in box(...): "plot.conf" n'est pas un paramètre graphique
#summary(mean)
checkresiduals(mean)
##
## Ljung-Box test
##
## data: Residuals from Mean
## Q* = 731.64, df = 18, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 19accuracy(mean, data_ts_test)# MAPE : 22%## ME RMSE MAE MPE MAPE MASE
## Training set 1.941958e-11 65611.93 55535.08 -3.855657 17.01186 2.161770
## Test set 1.037870e+05 110031.92 103787.04 22.486144 22.48614 4.040036
## ACF1 Theil's U
## Training set 0.8254447 NA
## Test set 0.0485288 2.517689summary(naivem)##
## Forecast method: Naive method
##
## Model Information:
## Call: naive(y = data_ts_train, h = 8)
##
## Residual sd: 36679.9508
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 1896.811 36679.95 29013.27 -0.02007386 8.597313 1.129377
## ACF1
## Training set -0.2744236
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2019 426097 379089.8 473104.2 354205.6 497988.4
## Feb 2019 426097 359618.7 492575.3 324427.2 527766.8
## Mar 2019 426097 344678.1 507515.9 301577.5 550616.5
## Apr 2019 426097 332082.5 520111.5 282314.2 569879.8
## May 2019 426097 320985.6 531208.4 265343.0 586851.0
## Jun 2019 426097 310953.2 541240.8 249999.8 602194.2
## Jul 2019 426097 301727.5 550466.5 235890.3 616303.7
## Aug 2019 426097 293140.4 559053.6 222757.5 629436.5checkresiduals(naivem)
##
## Ljung-Box test
##
## data: Residuals from Naive method
## Q* = 248.52, df = 19, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 19accuracy(naivem, data_ts_test) # MAPE : 8.5%## ME RMSE MAE MPE MAPE MASE
## Training set 1896.811 36679.95 29013.27 -0.02007386 8.597313 1.129377
## Test set 24357.125 43915.19 38328.62 4.72582164 8.499751 1.491988
## ACF1 Theil's U
## Training set -0.2744236 NA
## Test set 0.0485288 1.063155summary(driftm)##
## Forecast method: Random walk with drift
##
## Model Information:
## Call: rwf(y = data_ts_train, h = 8, drift = T)
##
## Drift: 1896.8105 (se 3778.1861)
## Residual sd: 36825.2032
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 2.297696e-11 36630.87 28899.04 -0.5861884 8.591266 1.12493
## ACF1
## Training set -0.2744236
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2019 427993.8 380800.4 475187.2 355817.7 500169.9
## Feb 2019 429890.6 362798.7 496982.5 327282.4 532498.8
## Mar 2019 431787.4 349190.1 514384.7 305465.7 558109.1
## Apr 2019 433684.2 337818.7 529549.8 287070.6 580297.9
## May 2019 435581.1 327854.7 543307.4 270827.8 600334.3
## Jun 2019 437477.9 318875.0 556080.7 256090.5 618865.2
## Jul 2019 439374.7 310630.0 568119.3 242476.7 636272.6
## Aug 2019 441271.5 302958.0 579585.0 229739.3 652803.7checkresiduals(driftm)
##
## Ljung-Box test
##
## data: Residuals from Random walk with drift
## Q* = 248.52, df = 18, p-value < 2.2e-16
##
## Model df: 1. Total lags used: 19accuracy(driftm, data_ts_test) # MAPE : 7.5%## ME RMSE MAE MPE MAPE MASE
## Training set 2.297696e-11 36630.87 28899.04 -0.5861884 8.591266 1.124930
## Test set 1.582148e+04 41314.60 33586.60 2.7843152 7.582963 1.307399
## ACF1 Theil's U
## Training set -0.27442358 NA
## Test set 0.06907259 1.007801summary(snaivem)##
## Forecast method: Seasonal naive method
##
## Model Information:
## Call: snaive(y = data_ts_train, h = 8)
##
## Residual sd: 28666.7301
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 25337.46 28666.73 25689.63 7.101745 7.207375 1 0.2695124
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2019 415292 378554.1 452029.9 359106.2 471477.8
## Feb 2019 423665 386927.1 460402.9 367479.2 479850.8
## Mar 2019 478207 441469.1 514944.9 422021.2 534392.8
## Apr 2019 443548 406810.1 480285.9 387362.2 499733.8
## May 2019 464162 427424.1 500899.9 407976.2 520347.8
## Jun 2019 457944 421206.1 494681.9 401758.2 514129.8
## Jul 2019 440436 403698.1 477173.9 384250.2 496621.8
## Aug 2019 366272 329534.1 403009.9 310086.2 422457.8checkresiduals(snaivem)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 35.426, df = 19, p-value = 0.0124
##
## Model df: 0. Total lags used: 19accuracy(snaivem, data_ts_test) # MAPE : 3.6%## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 25337.46 28666.73 25689.63 7.101745 7.207375 1.0000000 0.2695124
## Test set 14263.38 22148.43 16960.88 3.053421 3.648407 0.6602226 -0.5427745
## Theil's U
## Training set NA
## Test set 0.4792835Etude du Modèle de Buys-Ballot
Modèle
L’approche de BUYS-BALLOT consiste à introduire des variables indicatrices correspondant à chaque saison définit par le cycle d’observation. Pour les données trimestrielles, on intègre 4 variables indicatrices. Et pour les données mensuelles, on intègre 12 variables indicatrices.
Le modèle doit alors être estimé (sans constante) avec ces variables indicatrices.
Soit le modèle de Buys-Ballot tel que \(Xt = Zt + St + \mu t\)
Nous allons estimé dans un premier temp la tendance \(Zt\), puis dans un second temps la saisonnalité \(St\), tandis que les \(\mu t\) ne pourrons pas être estimer puisqu’il s’agit par définition d’accidents.
Préparation des données
Annees=as.numeric(time(data_ts_train))
ts_DataFrame =data.frame(trafic=data_ts_train,X=as.numeric(Annees))Estimation de la tendance
Regression <- lm(trafic~X,data = ts_DataFrame)On utilise la droite de régression obtenu pour prédire la tendance des 8 mois de 2019
tendance=predict(Regression)
#les 8 prochains mois
AnneeMoisNumericFutur=seq(max(Annees)+1/12,length=8,by=1/12)
tendance2=predict(Regression, newdata=data.frame(X=AnneeMoisNumericFutur)) Estimation de la saisonnalité
On récupère les résidus issu de l’estimation de la tendance.
ts_DataFrame$trafic_residual <- residuals(Regression)On définit la saisonnalité avec les mois.
ts_DataFrame$mois <- round(ts_DataFrame$X - trunc(ts_DataFrame$X),digit=4)On estime la saisonnalité à partir des résidus du modèle comportant uniquement la tendance.
Regression2 =lm(trafic_residual~0+as.factor(mois),data=ts_DataFrame)Prédiction de la saisonnalité sur le jeu d’entraînement
prediction2 =predict(Regression2)Prédiction sur les 8 mois de 2019
MoisNumeric= round(AnneeMoisNumericFutur - trunc(AnneeMoisNumericFutur
),4)
Prediction3 =predict( Regression2, newdata= data.frame(mois=MoisNumeric))Calculons une région de confiance avec l’erreur d’ajustement
ResidusRegression2=residuals(Regression2)
hist(ResidusRegression2)
1.96*sqrt(var(ResidusRegression2))## [1] 19226.16Auto corrélation de la série temporelle
L’auto corrélation de notre série temporelle correspond à la corrélation entre une mesure du trafic \(t\) et les mesures précédentes \(t - k\) ou les mesures suivantes \(t + k\).
L’auto covariance d’une variable \(Xt\) de moyenne \(\mu\) et d’écart type \(\sigma\) à un décalage \(k\) est donné par la formule
\(\gamma_k= E((X_t-\mu)(X_{t+k}-\mu))\)
On en déduit l’auto-corrélation correspondante :
\(\rho_k=\frac{\gamma_k}{\sigma^2}\)
Affichons les auto-corrélations de la séries prédite grâce à un corrélogramme
ACF_Sur_Valeurs_Predites <- acf(prediction2)
Il est normal que la série soit auto corrélé totalement à elle avec un décalage nulle.
On observe une corrélation forte (0.87) avec un décalage (lag) de 12, cela correspond bien à une saisonnalité annuelle.
print(data.frame(ACF_Sur_Valeurs_Predites$lag,ACF_Sur_Valeurs_Predites$acf)[1:13,])## ACF_Sur_Valeurs_Predites.lag ACF_Sur_Valeurs_Predites.acf
## 1 0 1.00000000
## 2 1 0.13556225
## 3 2 -0.30353869
## 4 3 0.19244891
## 5 4 -0.11749420
## 6 5 -0.36058135
## 7 6 -0.03501791
## 8 7 -0.32578902
## 9 8 -0.14275858
## 10 9 0.15718886
## 11 10 -0.25819898
## 12 11 0.12067802
## 13 12 0.87497657Recalculons la valeur d’auto-corrélation obtenu en appliquant la formule.
Observons l’application de la formule, en choisissant un décalage de 12
#Constantes
Nombre_Observations=96
decalage=12
#Estimations
moyenneMu=mean(prediction2)
sdSigma=sd(prediction2)
Serie1=prediction2[(decalage+1): 96 ]
Serie2=prediction2[ 1 :(96-decalage)]
GammaDecalage12=mean((Serie1-moyenneMu)*(Serie2-moyenneMu))*((Nombre_Observations-decalage)/(Nombre_Observations))
RhoDecalage12=GammaDecalage12/(sdSigma^2)
RhoDecalage12## [1] 0.8658622Le résultat obtenu est correct. L’auto corrélation avec un décalage de 12 est donc très forte.
De plus cette auto corrélation étant positive, cela indique une tendance croissante.
la deuxième plus forte corrélation est observé avec un décalage de 5, observons cela graphiquement
plot ( 1:length(prediction2), prediction2,type="l")
points((1:length(prediction2))-5,prediction2,type="l",col="red")
Cette corrélation est peu pertinente.
Après avoir étudier les auto-corrélations sur l’ensemble du modèle, Observons les auto-corrélations sur les résidus du modèle de Buys-Ballot.
Une auto corrélation de ces résidus signifierait que les accidents seraient prévisibles, cela n’aurait aucun sens et invaliderait notre modèle.
plot(acf(ResidusRegression2))
Pour notre modèle, il n’y a aucune auto-corrélation significative des résidus (symbolisé par la ligne bleu).
Notre modèle est donc viable.
Comparaison des prédictions et des valeurs réelles
Affichage de la tendance
Affichage du modèle de Buys Ballot
Affichage de la prédiction sur les 8 mois de 2020
Préparation DataFrame pour affichage ggplot
Reproduisons les graphiques avec ggplot2 pour un résultat plus professsionnel.
plotBuysBallot <- Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$BuysBalotModele)
plotBuysBallot
Sauvegarde de l’image pour utilisation ultérieur dans l’application web.
Nous avons réussi à Créer un modèle de prédiction correct avec Buys-Ballot. On remarque que la prédiction semble bien correspondre à la réalité si on fait abstraction du dernier mois où le nombre de voyageurs a bien plus chuté que la prédiction du modèle de Buys-Balot.
Comparons avec un ajustement local réalisé par lissage moyennes mobiles.
Comparaison avec les valeurs observées
Lissage moyenne mobile
Définition
Une moyenne mobile est un filtre linéaire. Il permet de transformer une série chronologique avec comme but d’annuler une composante (tendance ou saison) pour en laisser les autres invariantes tout en réduisant le bruit.
Une moyenne mobile en t est définit comme une combinaison linéaire finie des valeurs de la série correspondantes à des dates entourant t, c’est donc un lissage de la série.
Une moyenne mobile d’ordre m peut être écrite tel que \(\begin{equation} \hat{T}_{t} = \frac{1}{m} \sum_{j=-k}^k y_{t+j}, \end{equation}\)
Avec \(m=2k+1\)
L’estimation de la tendance + saisonnalité au temps t est obtenue en faisant la moyenne des valeurs de la série chronologique dans les k périodes de t. Les observations qui sont proches dans le temps sont également susceptibles d’être proches en valeur estimé.
Par conséquent, la moyenne élimine une partie du caractère aléatoire des données, laissant une composante tendance + saisonnalité lisse.
Choix Moyenne mobiles
Commençons notre analyse avec une moyenne mobile d’ordre 11, décidé de façon arbitraire
library(forecast)
MM_ordre_11 <- ma(DataAffichageGGplot$trafic[1:96] , 11,centre=FALSE)La première valeur (non manquante) de cette colonne est la moyenne des onze premières observations càd de janvier à novembre 2011. La deuxième valeur est la moyenne des valeurs de février 2011 à décembre 2011 et ainsi de suite. De façon à que chaque valeur est la moyenne des observations sur 11 mois centrée sur le mois correspondant.
Il est donc normale qu’il manque 5 estimations pour les premiers & derniers mois puisque les données antérieur à 2011 ne sont pas disponibles tout comme les données futur sont inconnus.
Il existe une méthode qui permet l’estimation des points extrêmes que nous utiliserons pour la suite.
Affichage de la moyenne mobile simple d’ordre 11
Affichage_Prediction(DataAffichageGGplot,c(MM_ordre_11, rep(NA, 8)) )## Warning: Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 104 row(s) containing missing values (geom_path).## Warning: Removed 18 row(s) containing missing values (geom_path).
On remarque que les données estimées (ici en bleu) représente bien la tendance générale de la série, dans le cas d’une moyenne mobile, plus son ordre est élevé, plus la courbe des estimations sera lisse.
Voici pour exemple l’affichage de différentes moyennes mobiles simples.
## Warning: Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 104 row(s) containing missing values (geom_path).## Warning: Removed 10 row(s) containing missing values (geom_path).
## Warning: Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 104 row(s) containing missing values (geom_path).## Warning: Removed 12 row(s) containing missing values (geom_path).
## Warning: Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 104 row(s) containing missing values (geom_path).## Warning: Removed 22 row(s) containing missing values (geom_path).
Nous avons choisi des moyennes mobiles simples d’un ordre impair pour qu’elles soient symétriques. En effet, dans une moyenne mobile d’ordre m=2k+1, on fait la moyenne de l’observation centrale et des k observations de chaque côté. Mais si m était pair, elle ne serait plus symétrique.
Moyenne mobile d’une Moyenne mobile
Il est possible d’appliquer une moyenne mobile à une moyenne mobile. Cela permet de rendre symétrique une moyenne mobile d’ordre pair.
Appliquons cela par une moyenne mobile d’ordre 2 sur une moyenne mobile d’ordre 4
MM_ordre_4 <- ma(DataAffichageGGplot$trafic[1:96] , 4,centre=FALSE)
MM_Ordre_2x4 <- ma(DataAffichageGGplot$trafic[1:96], 4, centre=TRUE)
head(MM_ordre_4)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] NA 256279.5 266477.2 282748.5 281141.8 267147.0head(MM_Ordre_2x4)## Time Series:
## Start = 1
## End = 6
## Frequency = 1
## [1] NA NA 261378.4 274612.9 281945.1 274144.4Une MM 2x4 signifie une MM d’ordre 4 suivi d’une MM d’ordre 2. Concrètement, les 5 premières valeurs de la série chronologique sont 245900, 238000, 263227,277991 et 286691. Les deux premiers termes de la MM d’ordre 4 sera
(245900+238000+263227+277991) /4 = 1 025 118/4 = 256279.5
(238000+263227+277991+286691) /4 = 266477.2
Donc le premier terme de la moyenne mobile d’ordre 2, calculé à partir des 2 premiers termes de la MM d’ordre 4, donnera (256279.5 + 266477.2)/2 = 261378.4
On combinant ces deux moyennes mobiles, nous avons obtenu une moyenne mobile centrée d’ordre 4, les valeurs calculées obtenues sont symétriques, en effet, nous avons appliqué cela aux données :
\(\begin{align*} \hat{T}_{t} &= \frac{1}{2}\Big[ \frac{1}{4} (y_{t-2}+y_{t-1}+y_{t}+y_{t+1}) + \frac{1}{4} (y_{t-1}+y_{t}+y_{t+1}+y_{t+2})\Big] \\ &= \frac{1}{8}y_{t-2}+\frac14y_{t-1} + \frac14y_{t}+\frac14y_{t+1}+\frac18y_{t+2}. \end{align*}\)
Ainsi si l’on combine 2 moyennes mobiles, nous devons appliqués une moyenne mobile pair après une moyenne mobile pair et une moyenne mobile impair après une moyenne mobile impair afin que nos résultats soit symétriques.
Conservation & Annulation
Soit la Moyenne mobile 2x12 \(\hat{T}_{t} = \frac{1}{32}y_{t-6} + \frac{1}{12}y_{t-5} + \frac{1}{12}y_{t-4} + \frac{1}{12}y_{t-3} + \frac{1}{12}y_{t-2} + \frac{1}{12}y_{t-1} +\frac{1}{12}y_{t} + \frac{1}{12}y_{t+1} +\frac{1}{12}y_{t+2} + \frac{1}{12}y_{t+3} + \frac{1}{12}y_{t+4} + \frac{1}{12}y_{t+5} + \frac{1}{32}y_{t+6}.\)
Appliquée au trafic qui a une saisonnalité annuelle, chaque mois de l’année a le même poids, car le premier et le dernier terme s’appliquent au même mois des années consécutives. Ainsi, la variation saisonnière est éliminée.
Un autre choix de Moyenne mobile entraînera une estimation soumis à la saisonnalité.
MM_ordre_2x12 <- ma(DataAffichageGGplot$trafic[1:96], 12)
Affichage_Prediction(DataAffichageGGplot,c(MM_ordre_2x12, rep(NA, 8)) ) + labs(title = "Moyenne Mobile d'ordre 2x12")## Warning: Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 104 row(s) containing missing values (geom_path).## Warning: Removed 20 row(s) containing missing values (geom_path).
La droite de prédiction obtenu ne présente aucune trace de saisonnalité. Un autre choix de moyenne mobile (hormis 12k \(k \in N\)) aurait présenter des traces de saisonnalités.
Nous avons vu que la combinaison de moyennes mobiles produisait des moyennes mobiles pondérés, la combinaison que nous avons utilisés produisait les poids suivants \(\left[\frac{1}{32},\frac{1}{32},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{12},\frac{1}{32}\right]\)
En effet, une moyenne mobile d’ordre m s’écrit tel que \(\hat{T}_t = \sum_{j=-k}^k a_j y_{t+j}\) Avec \(k=(m-1)/2\) et les poids \(\left[a_{-k},\dots,a_k\right]\) Si leur somme vaut 1 et que ces poids sont symétriques (tq \(a_j = a_{-j}\)) alors la courbe obtenu sera plus lisse.
En conclusion, nous avons pu éliminer les saisonnalités de période 12 et conserver la tendance, notre choix d’ordre de moyenne mobile est correct, on remarque que choisir un ordre arbitrairement grand est inutile, il suffit de comparer les moyennes mobiles d’ordre 15 et 12, celle à 12 annule de façon bien plus performante les saisonnalités.
Lissage exponentielle
L’idée de la méthode de lissage exponentielle est une extension de la méthode naîve (naivm vu précédemment). Le principe étant de fournir une prévision qui dépend de moyenne pondérées d’observations passées, les poids diminuant de façon exponentielle à mesure que les observations vieillissent.
Plus simplement, les observations les plus récentes ont un poids plus important que les observations les plus anciennes.
Le paramètre de lissage est décidé par le paramètre α.
On remarque que pour un lissage simple, la prédiction est égale à la valeur de la dernière observation.
Lissage simple
fcst_se <- ses(data_ts_train, h = 8)
summary(fcst_se)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = data_ts_train, h = 8)
##
## Smoothing parameters:
## alpha = 0.2559
##
## Initial states:
## l = 258126.0245
##
## sigma: 31480.96
##
## AIC AICc BIC
## 2430.727 2430.988 2438.420
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 7057.127 31151.31 25752.89 1.326234 7.684321 1.002462 0.0711143
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2019 431488.3 391143.9 471832.8 369786.8 493189.9
## Feb 2019 431488.3 389843.9 473132.7 367798.7 495178.0
## Mar 2019 431488.3 388583.3 474393.3 365870.8 497105.9
## Apr 2019 431488.3 387358.8 475617.9 363998.0 498978.7
## May 2019 431488.3 386167.3 476809.4 362175.7 500800.9
## Jun 2019 431488.3 385006.3 477970.4 360400.2 502576.5
## Jul 2019 431488.3 383873.6 479103.1 358667.9 504308.8
## Aug 2019 431488.3 382767.3 480209.4 356975.9 506000.8checkresiduals(fcst_se)
##
## Ljung-Box test
##
## data: Residuals from Simple exponential smoothing
## Q* = 144.66, df = 17, p-value < 2.2e-16
##
## Model df: 2. Total lags used: 19plot(fcst_se)
lines(data_ts_test, col="red")
df_se = as.data.frame(fcst_se)
predict_value_se <- df_se$`Point Forecast`
MAPE(predict_value_se, data_ts_test) # MAPE = 7.95## [1] 7.945782Optimisation du modèle
La librairie Forecast permet de trouver automatiquement le meilleur degré de lissage exponentielle fournissant la meilleure prédiction.
fit_ets <- ets(data_ts_train)
print(summary(fit_ets))## ETS(A,A,A)
##
## Call:
## ets(y = data_ts_train)
##
## Smoothing parameters:
## alpha = 0.1568
## beta = 1e-04
## gamma = 1e-04
##
## Initial states:
## l = 248267.1099
## b = 2163.3982
## s = -17928.3 -29535.73 9295.935 11005.81 -57117.85 -7708.17
## 38272.64 14592.34 16899.53 34763.15 -7344.204 -5195.15
##
## sigma: 11014.45
##
## AIC AICc BIC
## 2241.611 2249.458 2285.205
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -458.6799 10054.77 7831.554 -0.2623253 2.371375 0.3048527
## ACF1
## Training set 0.09626331checkresiduals(fit_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,A,A)
## Q* = 7.1794, df = 3, p-value = 0.06639
##
## Model df: 16. Total lags used: 19fcst_ets <- forecast(fit_ets, h=8)
plot(fcst_ets)
lines(data_ts_test, col="red")
df_ets = as.data.frame(fcst_ets)
predict_value_ets = df_ets$`Point Forecast`
MAPE(predict_value_ets, data_ts_test) #MAPE = 2.84## [1] 2.839718Affichage GGplot
DataAffichageGGplot$ModeleLissageExponentielle <- c(fcst_ets$fitted ,predict_value_ets )
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ModeleLissageExponentielle)## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).
Modèle ARIMA / SAMIRA Automatique
ARIMA : AutoRegressive Integrated Moving Average
Le modèle ARIMA est une combinaison du modèle ARMA* combiné à une différentiation (le Integrated)
*Le modèle ARMA (Autoregressive Moving Average Model) est définit par l’idée de décrire une séries temporelles en deux polynômes. Le premier étant pour l’auto regréssion et le second pour la moyenne mobile. ARMA(p,q) où :
p est l’ordre du polynôme auto régressif,
q est l’ordre du polynôme de moyenne mobile L’équation est données par avec :
φ = les paramètres du modèle autorégressif,
θ = les paramètres du modèle de moyenne mobile.
c = une constante,
Σ = notation de sommation,
ε = termes d’erreur (bruit blanc).
Une différentiation correspond à retirer les tendances
tendance linéaire : une différenciation
tendance quadratique : deux différenciations
Le modèle SARIMA est une combinaison du modèle ARIMA qui prend en compte la composante saisonnière.
Auto.arima prend en compte les saisonnalités, comme on peut le voir dans le modèle sélectionné : (0,1,1)(0,1,1)[12]
# retourne les meilleurs paramètres
# d=1 enleve la tendance
# D=1 enleve la saisonnalité
# => avoir des données stationnaires
# trace : voir les résultats
fit_arima <- auto.arima(data_ts_train, d=1, D=1, stepwise = FALSE, approximation = FALSE, trace=TRUE)##
## ARIMA(0,1,0)(0,1,0)[12] : 1846.398
## ARIMA(0,1,0)(0,1,1)[12] : 1833.134
## ARIMA(0,1,0)(0,1,2)[12] : 1835.211
## ARIMA(0,1,0)(1,1,0)[12] : 1833.056
## ARIMA(0,1,0)(1,1,1)[12] : 1835.09
## ARIMA(0,1,0)(1,1,2)[12] : Inf
## ARIMA(0,1,0)(2,1,0)[12] : 1835.207
## ARIMA(0,1,0)(2,1,1)[12] : 1837.012
## ARIMA(0,1,0)(2,1,2)[12] : 1836.461
## ARIMA(0,1,1)(0,1,0)[12] : 1814.951
## ARIMA(0,1,1)(0,1,1)[12] : 1801.155
## ARIMA(0,1,1)(0,1,2)[12] : 1803.362
## ARIMA(0,1,1)(1,1,0)[12] : 1803.592
## ARIMA(0,1,1)(1,1,1)[12] : 1803.361
## ARIMA(0,1,1)(1,1,2)[12] : Inf
## ARIMA(0,1,1)(2,1,0)[12] : 1805.004
## ARIMA(0,1,1)(2,1,1)[12] : 1805.397
## ARIMA(0,1,1)(2,1,2)[12] : Inf
## ARIMA(0,1,2)(0,1,0)[12] : 1816.915
## ARIMA(0,1,2)(0,1,1)[12] : 1803.033
## ARIMA(0,1,2)(0,1,2)[12] : 1805.296
## ARIMA(0,1,2)(1,1,0)[12] : 1805.702
## ARIMA(0,1,2)(1,1,1)[12] : 1805.295
## ARIMA(0,1,2)(1,1,2)[12] : Inf
## ARIMA(0,1,2)(2,1,0)[12] : 1807.026
## ARIMA(0,1,2)(2,1,1)[12] : 1807.441
## ARIMA(0,1,3)(0,1,0)[12] : 1817.787
## ARIMA(0,1,3)(0,1,1)[12] : Inf
## ARIMA(0,1,3)(0,1,2)[12] : Inf
## ARIMA(0,1,3)(1,1,0)[12] : Inf
## ARIMA(0,1,3)(1,1,1)[12] : Inf
## ARIMA(0,1,3)(2,1,0)[12] : Inf
## ARIMA(0,1,4)(0,1,0)[12] : 1820.052
## ARIMA(0,1,4)(0,1,1)[12] : Inf
## ARIMA(0,1,4)(1,1,0)[12] : Inf
## ARIMA(0,1,5)(0,1,0)[12] : Inf
## ARIMA(1,1,0)(0,1,0)[12] : 1825.579
## ARIMA(1,1,0)(0,1,1)[12] : 1812.512
## ARIMA(1,1,0)(0,1,2)[12] : 1814.657
## ARIMA(1,1,0)(1,1,0)[12] : 1813.2
## ARIMA(1,1,0)(1,1,1)[12] : 1814.614
## ARIMA(1,1,0)(1,1,2)[12] : 1816.192
## ARIMA(1,1,0)(2,1,0)[12] : 1815.227
## ARIMA(1,1,0)(2,1,1)[12] : Inf
## ARIMA(1,1,0)(2,1,2)[12] : 1817.796
## ARIMA(1,1,1)(0,1,0)[12] : 1816.841
## ARIMA(1,1,1)(0,1,1)[12] : 1802.853
## ARIMA(1,1,1)(0,1,2)[12] : 1805.117
## ARIMA(1,1,1)(1,1,0)[12] : 1805.653
## ARIMA(1,1,1)(1,1,1)[12] : Inf
## ARIMA(1,1,1)(1,1,2)[12] : Inf
## ARIMA(1,1,1)(2,1,0)[12] : Inf
## ARIMA(1,1,1)(2,1,1)[12] : Inf
## ARIMA(1,1,2)(0,1,0)[12] : 1819.234
## ARIMA(1,1,2)(0,1,1)[12] : 1805.22
## ARIMA(1,1,2)(0,1,2)[12] : 1807.539
## ARIMA(1,1,2)(1,1,0)[12] : 1807.381
## ARIMA(1,1,2)(1,1,1)[12] : 1807.538
## ARIMA(1,1,2)(2,1,0)[12] : 1808.925
## ARIMA(1,1,3)(0,1,0)[12] : 1820.05
## ARIMA(1,1,3)(0,1,1)[12] : 1806.055
## ARIMA(1,1,3)(1,1,0)[12] : 1808.732
## ARIMA(1,1,4)(0,1,0)[12] : Inf
## ARIMA(2,1,0)(0,1,0)[12] : 1824.435
## ARIMA(2,1,0)(0,1,1)[12] : 1811.07
## ARIMA(2,1,0)(0,1,2)[12] : 1813.287
## ARIMA(2,1,0)(1,1,0)[12] : 1811.619
## ARIMA(2,1,0)(1,1,1)[12] : 1813.247
## ARIMA(2,1,0)(1,1,2)[12] : Inf
## ARIMA(2,1,0)(2,1,0)[12] : 1813.821
## ARIMA(2,1,0)(2,1,1)[12] : 1815.872
## ARIMA(2,1,1)(0,1,0)[12] : Inf
## ARIMA(2,1,1)(0,1,1)[12] : Inf
## ARIMA(2,1,1)(0,1,2)[12] : Inf
## ARIMA(2,1,1)(1,1,0)[12] : Inf
## ARIMA(2,1,1)(1,1,1)[12] : Inf
## ARIMA(2,1,1)(2,1,0)[12] : Inf
## ARIMA(2,1,2)(0,1,0)[12] : Inf
## ARIMA(2,1,2)(0,1,1)[12] : Inf
## ARIMA(2,1,2)(1,1,0)[12] : Inf
## ARIMA(2,1,3)(0,1,0)[12] : Inf
## ARIMA(3,1,0)(0,1,0)[12] : 1823.646
## ARIMA(3,1,0)(0,1,1)[12] : 1808.49
## ARIMA(3,1,0)(0,1,2)[12] : 1810.542
## ARIMA(3,1,0)(1,1,0)[12] : 1808.594
## ARIMA(3,1,0)(1,1,1)[12] : 1810.321
## ARIMA(3,1,0)(2,1,0)[12] : 1810.708
## ARIMA(3,1,1)(0,1,0)[12] : Inf
## ARIMA(3,1,1)(0,1,1)[12] : Inf
## ARIMA(3,1,1)(1,1,0)[12] : Inf
## ARIMA(3,1,2)(0,1,0)[12] : Inf
## ARIMA(4,1,0)(0,1,0)[12] : 1823.996
## ARIMA(4,1,0)(0,1,1)[12] : 1810.199
## ARIMA(4,1,0)(1,1,0)[12] : 1810.845
## ARIMA(4,1,1)(0,1,0)[12] : Inf
## ARIMA(5,1,0)(0,1,0)[12] : 1825.055
##
##
##
## Best model: ARIMA(0,1,1)(0,1,1)[12]print(summary(fit_arima))## Series: data_ts_train
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.7675 -0.5465
## s.e. 0.0977 0.1295
##
## sigma^2 = 138827719: log likelihood = -897.43
## AIC=1800.85 AICc=1801.16 BIC=1808.11
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 805.2378 10822.93 7747.841 0.2158443 2.213481 0.3015941 0.03453594checkresiduals(fit_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 10.898, df = 17, p-value = 0.8618
##
## Model df: 2. Total lags used: 19fcst_arima <- forecast(fit_arima, h=8)
plot(fcst_arima)
lines(data_ts_test, col='red')
df_arima = as.data.frame(fcst_arima)
predict_value_arima = df_arima$`Point Forecast`
MAPE(predict_value_arima, data_ts_test)## [1] 2.72507Affichage GGplot
DataAffichageGGplot$ModeleArima <- c(fit_arima$fitted ,predict_value_arima )
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ModeleArima)## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).
Prophet
Prophet est un outil de prévision de série chronologique open-source mis en production par Facebook.
L’idée est décomposer la série temporelle sous la forme :
y(t) = g(t) + s(t) + h(t) + e(t)
avec respectivement :
- g(t) : la tendance (linéaire ou logistique)
- s(s) : une ou plusieurs composantes saisonnières (annuelle, hebdomadaire ou quotidienne)
- h(t) : l’effet des vacances ou de jours spécifiques qui pourront être paramétrés
- e(t) : l’erreur, bruit aléatoire qui mesure l’écart entre le modèle et les données réelles
Ce modèle est particulièrement utile quand on a de forte saisonnalités et permet de prendre en compte les vacances (et ses effets) par pays.
Préparation des données
format des dates
colonne du temps name : ds
colonne de la mesure quantitative name : y
library(prophet)## Le chargement a nécessité le package : Rcpp## Le chargement a nécessité le package : rlang##
## Attachement du package : 'rlang'## Les objets suivants sont masqués depuis 'package:purrr':
##
## %@%, as_function, flatten, flatten_chr, flatten_dbl, flatten_int,
## flatten_lgl, flatten_raw, invoke, splicelibrary(zoo)##
## Attachement du package : 'zoo'## Les objets suivants sont masqués depuis 'package:base':
##
## as.Date, as.Date.numericdata_train$ds <- as.Date( as.yearmon(time(data_ts_train)))On créer notre modèle et on prédit sur les 8 mois de 2019.
On affiche les courbes.
On regarde la valeur de la MAPE : 3.2.
model_prophet <- prophet(data_train)## Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.## Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.forecast_prophet <- make_future_dataframe(model_prophet, periods = 8, freq = 'month')
AAPLfc <- predict(model_prophet, forecast_prophet)
tail(AAPLfc[c("ds", "yhat", "yhat_lower", "yhat")])## ds yhat yhat_lower yhat.1
## 99 2019-03-01 493929.1 482512.9 493929.1
## 100 2019-04-01 477074.4 465193.4 477074.4
## 101 2019-05-01 475880.5 462887.3 475880.5
## 102 2019-06-01 503278.4 491122.3 503278.4
## 103 2019-07-01 457437.7 445628.0 457437.7
## 104 2019-08-01 412673.2 401309.4 412673.2dyplot.prophet(model_prophet, AAPLfc)## Warning: `select_()` was deprecated in dplyr 0.7.0.
## Please use `select()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.data_pp <- subset(AAPLfc, select=c("yhat"))
data_pp_ts <- ts(data_pp, start=2011, frequency=12)
data_pp_ts_w <- window(data_pp_ts, start= c(2019,1), end = c(2019,8))
MAPE(data_pp_ts_w, data_ts_test) #3.2## [1] 3.189021Affichage GGplot
DataAffichageGGplot$ProphetModele <- AAPLfc$yhat
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$ProphetModele)## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).
LSTM
Un réseau LSTM (Long Short Term Memory) est un type particulier de réseau de neurones récurrent qui est capable d’apprendre les dépendances à long terme. Il a été introduit pas Hochreiter & Schmidhuber en 1997.
L’idée d’un réseaux récurrent est de mettre en place des boucles qui permettent aux informations de persister. Ils sont représentés par une chaîne de modules répétitifs de réseau de neurones.
Chaque module est traversé un à un en faisant évoluer son état qui peut se voir ajouter ou supprimer de l’information, ces changements sont régulés par ce qu’on appelle des “portes”.
Elles sont définit par une couche de neurones sigmoïde et d’opération de multiplication ponctuelle qui laisse passer ou non l’information. La fonction Sigmoïde produit des nombres entre 1 et 0 ce qui définit la quantité de chaque composant qui doit passer.
Préparation des données
On veut prendre l’année précédente pour apprendre , on fixe donc lag à 12, en réalité cela correspond à faire 12 - 1.
Puis on transforme en array 3D car le modèle LSTM prend un tensor de format 3D [samples, timesteps, features]
samples : nbr d’observation par batchs
timesteps : lag
features : nbr de valeur predites
Construction de notre réseaux de neuronnes
lstm_model <- keras_model_sequential()## Loaded Tensorflow version 2.8.0lstm_model %>%
layer_lstm(units = 50, # taille du layer
batch_input_shape = c(1, 12, 1), # batch size, timesteps, features
return_sequences = TRUE,
stateful = TRUE) %>%
# pourcentage de neuronnes qu'on détruit à chaque itération
layer_dropout(rate = 0.5) %>%
layer_lstm(units = 50,
return_sequences = TRUE,
stateful = TRUE) %>%
layer_dropout(rate = 0.5) %>%
time_distributed(keras::layer_dense(units = 1))
lstm_model %>%
compile(loss = 'mae', optimizer = 'adam', metrics = 'accuracy')
summary(lstm_model)## Model: "sequential"
## ________________________________________________________________________________
## Layer (type) Output Shape Param #
## ================================================================================
## lstm_1 (LSTM) (1, 12, 50) 10400
## dropout_1 (Dropout) (1, 12, 50) 0
## lstm (LSTM) (1, 12, 50) 20200
## dropout (Dropout) (1, 12, 50) 0
## time_distributed (TimeDistributed) (1, 12, 1) 51
## ================================================================================
## Total params: 30,651
## Trainable params: 30,651
## Non-trainable params: 0
## ________________________________________________________________________________Apprentissage sur le jeu d’entraînement
lstm_model %>% fit(
x = x_train_arr,
y = y_train_arr,
batch_size = 1,
epochs = 20,
verbose = 0,
shuffle = FALSE #pour preserver nos séquences
)Prédiction sur le jeu de Test
On a 50 résultats / prédictions par input donc on transforme pour une seule prédiction
lstm_forecast <- lstm_model %>%
predict(x_pred_arr, batch_size = 1) %>%
.[, , 1]
# rescale en format basique
lstm_forecast <- lstm_forecast * scale_factors[2] + scale_factors[1]
lstm_forecast
fitted <- predict(lstm_model, x_train_arr, batch_size = 1) %>%
.[, , 1]
if (dim(fitted)[2] > 1) {
fit <- c(fitted[, 1], fitted[dim(fitted)[1], 2:dim(fitted)[2]])
} else {
fit <- fitted[, 1]
}
# rescale final de nos données
fitted <- fit * scale_factors[2] + scale_factors[1]
fitted
fitted <- c(rep(NA, lag), fitted)
fitted
length(fitted)Affichage et évaluation
Nos résultats ne sont pas très convaincant même en faisant varier nos paramètres, ceci peut être explicable par le manque d’observations.
## Jan Feb Mar Apr May Jun Jul Aug
## 2019 469206.2 471857.1 437201.1 446503.8 458682.0 457769.1 478484.8 472561.8## Jan Feb Mar Apr May Jun Jul Aug
## 2019 443700 441499 480649 463680 453372 505190 445332 370211
Affichage GGplot
DataAffichageGGplot$LSTM_Modele <- fitted
Affichage_Prediction(DataAffichageGGplot, DataAffichageGGplot$LSTM_Modele)## Warning: Removed 98 row(s) containing missing values (geom_path).
## Removed 98 row(s) containing missing values (geom_path).## Warning: Removed 12 row(s) containing missing values (geom_path).
Comparaison Modèle
Choix des indicateurs
Nous devons pouvoir comparer nos modèles entre eux, nous souhaitons pouvoir observer leurs performances globale, ainsi que spécifiquement leurs performances sur le jeu d’entraînement et le jeu de test (les 8 mois de 2019 à prédire).
Nous utiliserons une mesure largement utilisé lors de nos précédents travaux : le \(R^2 \in [0 \; ; 1]\) , cette mesure indique la proportion de la variance expliquée par le modèle.
0 % le modèle n’explique par la variable Y
100 % le modèle explique la variabilité de Y lié à la liaison linéaire des variables explicatives entièrement
Rappelons son explication mathématiques , soit SCE la somme des distances au carré entre chaque valeur prédite par le modèle \({\widehat y_i}\) et la moyenne des réponses \(\overline{y}\)
\(\text{SCE} = \sum_{i=1}^{N}(\hat{y_i} – \overline{y})^2\)
Nous obtenons alors la part de dispersion expliquée par le modèle.
Puis, nous calculons la dispersion totale des données nommé SCT
\(\text{SCT} = \sum_{i=1}^{N}(y_{i } – \overline{y})^2\)
Avec \(y_i\) une valeur prise par une variable expliquée
Nous obtenons alors le R² par la combinaisons des calculs précédents
\(R^2 = \frac{SCE}{SCT}\)
Pour compléter cette mesure, nous utiliserons donc l’erreur absolue moyenne en pourcentage (MAPE en anglais) Il s’agit de la moyenne des écarts en valeur absolue par rapport aux valeurs observées.
C’est donc un pourcentage et par conséquent un indicateur pratique de comparaison.
Exemple de mise en production
Imaginons que l’on utilise nos modèles pour prédire le nombre de voyageurs, si l’on prévoit trop de passagers alors les moyens mis en place pour les accueillir sont partiellement utilisés, il en résulte un coût de 5 euros par passagers.
Si l’on a prédit moins de passagers que la réalité, notre compagnie doit commander en livraison urgente, il en résulte un coût de 10 euros par passagers.
Observons cette mise en production sur nos modèles.
Calculs des performances
Dans les concerts, chaque spectateur chante comme une casserole, pourtant la foule chante systématiquement juste, suivant cet intuition réalisons une prédiction qui est une moyenne des différents modèles réalisés.
DataAffichageGGplot$MoyenneDesModele <- rowMeans(DataAffichageGGplot[,5:9])Performances Modèles
noms_modeles <- c("Buys Ballot ","Lissage exponentielle ","Arima","Prophet","LSTM","Moyenne Des Modèles")
for (x in seq_along(noms_modeles)){
print("----------------------------")
print(noms_modeles[x])
predictionglobale <- DataAffichageGGplot[13:104,4+x]
predictionEntrainement <- DataAffichageGGplot[13:96,4+x]
predictionTest <- DataAffichageGGplot[97:104,4+x]
cat("MAPE globale : ",MAPE(DataAffichageGGplot$trafic[13:104], predictionglobale),"\n")
cat("R carré globale : ",Rcarre(DataAffichageGGplot$trafic[13:104], predictionglobale),"\n")
cat("MAPE Entrainement : ",MAPE(DataAffichageGGplot$trafic[13:96], predictionEntrainement),"\n")
cat("R carré Entrainement : ",Rcarre(DataAffichageGGplot$trafic[13:96], predictionEntrainement),"\n")
cat("MAPE Test : ",MAPE(DataAffichageGGplot$trafic[97:104], predictionTest),"\n")
cat("R carré Test : ",Rcarre(DataAffichageGGplot$trafic[97:104], predictionTest),"\n")
MiseEnProductionDuModele <- CoutDesErreurs(DataAffichageGGplot$trafic[97:104], predictionTest)
cat("Nombre de Passagers prévus en plus : ",MiseEnProductionDuModele[1],"\nNombre de Passagers prévus en moins : ",MiseEnProductionDuModele[2],"\nCout des erreurs en Euros : ",MiseEnProductionDuModele[3] )
cat("\nDifférence en nbre de voyageurs sur le dernier mois : ",predictionglobale[92] - DataAffichageGGplot$trafic[104],"\n")
}## [1] "----------------------------"
## [1] "Buys Ballot "
## MAPE globale : 2.205698
## R carré globale : 0.9740106
## MAPE Entrainement : 2.142334
## R carré Entrainement : 0.9747808
## MAPE Test : 2.871024
## R carré Test : 0.8036533
## Nombre de Passagers prévus en plus : 4848.58
## Nombre de Passagers prévus en moins : -91908.62
## Cout des erreurs en Euros : 1016058
## Différence en nbre de voyageurs sur le dernier mois : 36037.92
## [1] "----------------------------"
## [1] "Lissage exponentielle "
## MAPE globale : 2.270427
## R carré globale : 0.9725982
## MAPE Entrainement : 2.200386
## R carré Entrainement : 0.9736448
## MAPE Test : 3.005848
## R carré Test : 0.7863116
## Nombre de Passagers prévus en plus : 5276.056
## Nombre de Passagers prévus en moins : -96021.53
## Cout des erreurs en Euros : 1065736
## Différence en nbre de voyageurs sur le dernier mois : 38627.88
## [1] "----------------------------"
## [1] "Arima"
## MAPE globale : 2.551599
## R carré globale : 0.9655633
## MAPE Entrainement : 2.526595
## R carré Entrainement : 0.9629896
## MAPE Test : 2.814135
## R carré Test : 0.8420685
## Nombre de Passagers prévus en plus : 17397.14
## Nombre de Passagers prévus en moins : -81186.49
## Cout des erreurs en Euros : 1159808
## Différence en nbre de voyageurs sur le dernier mois : 24913.57
## [1] "----------------------------"
## [1] "Prophet"
## MAPE globale : 2.0477
## R carré globale : 0.9744149
## MAPE Entrainement : 1.919235
## R carré Entrainement : 0.9777721
## MAPE Test : 3.396584
## R carré Test : 0.7327864
## Nombre de Passagers prévus en plus : 1911.595
## Nombre de Passagers prévus en moins : -112557.7
## Cout des erreurs en Euros : 1163809
## Différence en nbre de voyageurs sur le dernier mois : 42462.18
## [1] "----------------------------"
## [1] "LSTM"
## MAPE globale : 3.439605
## R carré globale : 0.944196
## MAPE Entrainement : 3.5023
## R carré Entrainement : 0.9383416
## MAPE Test : 2.781304
## R carré Test : 0.7919626
## Nombre de Passagers prévus en plus : 48365.29
## Nombre de Passagers prévus en moins : -51742.04
## Cout des erreurs en Euros : 1484726
## Différence en nbre de voyageurs sur le dernier mois : 19874.38
## [1] "----------------------------"
## [1] "Moyenne Des Modèles"
## MAPE globale : 2.210271
## R carré globale : 0.9741124
## MAPE Entrainement : 2.151987
## R carré Entrainement : 0.9742663
## MAPE Test : 2.822246
## R carré Test : 0.8218652
## Nombre de Passagers prévus en plus : 12776.57
## Nombre de Passagers prévus en moins : -83900.11
## Cout des erreurs en Euros : 1094533
## Différence en nbre de voyageurs sur le dernier mois : 32383.19Sauvegardons le modèle de Buys Ballot pour une utilisation ultérieur dans une application Web
saveRDS(Regression, file = "./Shiny/RegressionAnnees.rda")
saveRDS(Regression2, file = "./Shiny/RegressionMois.rda")